Code
sshfs narecco@ant-login.linux.crg.es:/users/mirimia ~/mntFigure 4 and supplementary figure 4
Last code execution: 2024 February 07, Wednesday @ 19:02:17.
1 Intro
We performed histone PTMs DIA mass spectrometry and I analysed the raw files with EpiProfile v2.1 (Yuan et al. 2018). Here I parse the EpiProfile output and analyse the PTMs values with DEP (Zhang et al. 2018).
Acknowledgements: I’d like to thank Zuo-Fei Yuan and Simone Sidoli from the Ben Garcia lab for their help with data analysis and sample preparation respectively. Also a big thank you to Cristina Chiva from the CRG/UPF proteomics facility.
2 Set Up
To fetch files on the CRG cluster I mount the server on my local computer with the following command in the terminal.
2.1 Packages
Load packages required for the analysis and suppress any message. Check the Section 6 section at the end for package versions.
Code
library(dplyr, warn.conflicts = F, quietly = T)
library(readr)
library(readxl)
library(pheatmap)
library(DEP)
library(Cairo)
library(DT)Load my R package niar to use my custom made function to parse EpiProfile output.
Code
if ( ('niar' %in% .packages(all.available = TRUE)) == TRUE ) {
library('niar')
} else if ( ('niar' %in% .packages(all.available = TRUE)) == FALSE ){
message("The package niar is not installed so I'm trying to load it from ",
"the local repository of the package")
if ( dir.exists('~/mnt/narecco/software/R/niar' ) ) {
devtools::load_all(path = '~/mnt/narecco/software/R/niar')
} else {
stop("Can't find the local repo of the niar package! ",
"You must install it with:\n",
"devtools::install_github('Ni-Ar/niar') ")
}
} else{
stop("Can't understand if 'niar' package was installed beforehand")
}The package niar is not installed so I'm trying to load it from the local repository of the packageℹ Loading niar2.2 Functions
Define a plot style.
Code
theme_classic(base_size = 6.3, base_family = "Arial") +
theme(plot.title = element_text(size = 8, vjust = -1, hjust = 0),
plot.background = element_blank(),
strip.background = element_blank(),
panel.background = element_blank(),
axis.line = element_line(linewidth = 0.1),
axis.title.x = element_blank(),
axis.title.y = element_text(size = 5),
axis.text = element_text(colour = "black"),
axis.text.x = element_text(margin = margin(t = 0, unit = "mm")),
axis.ticks.x = element_blank(),
axis.ticks.y = element_line(linewidth = 0.1),
axis.ticks.length.y = unit(1, "mm"),
legend.margin = margin(l = 0, t = 0, unit = "mm"),
legend.box.margin = margin(l = 0, unit = "mm"),
legend.box.background = element_blank(),
legend.box.spacing = unit(0.5, "mm"),
# legend.title = element_blank(),
legend.text = element_text(margin = margin(l = -10, unit = "mm")),
legend.text.align = 0,
legend.key.height = unit(2, "mm"),
legend.key.width = unit(7.5, "mm"),
legend.key = element_blank(),
legend.background = element_blank(),
legend.spacing.x = unit(3, 'mm')
) -> alluvial_thm I modify the geom_stratum() function just to add a very small line width the the black borders of the columns. With linewidth = .05, I don’t have to change them manually later in Illustrator. This will hopefully fixed in future
Small function to parse the PTM abundance as 10 to the power scientific notion.
Code
scientific_10 <- function(x) {
parse(text=gsub("e\\+", " %*%10^", scales::scientific_format()(x)))
}2.3 Directories & File Paths
Here I organise all the file and directories paths I need to run the analysis and define where to save the processed tables and figures.
Code
oneDrive_Dir <- file.path("~/OneDrive - CRG - Centre de Regulacio Genomica/Suz12_AS_project")
code_dir_fig <- file.path(oneDrive_Dir, "_Code/Fig4")
tbl_dir_fig4 <- file.path(code_dir_fig, "tables")
pdf_dir_fig4 <- file.path(code_dir_fig, "pdfs")
if (!dir.exists(pdf_dir_fig4)) { dir.create(pdf_dir_fig4, recursive = T) }
if (!dir.exists(tbl_dir_fig4)) { dir.create(tbl_dir_fig4, recursive = T) }Path to medata data and EpiProfile output files.
Code
metadata_path <- file.path(tbl_dir_fig4, 'Histone_PTMs_metadata.xlsx')
WT_Dex4_histone_ratio_path <- list.files(tbl_dir_fig4,
pattern = "FullSet1_Run2_DIA_histone_ratios_Suz12dExon4.xls",
full.names = T)
Resce_KO_histone_ratio_path <- list.files(tbl_dir_fig4,
pattern = "FullSet2_Run1_DIA_histone_ratios_Suz12KO.xls",
full.names = T)
stopifnot(file.exists(metadata_path))
stopifnot(file.exists(WT_Dex4_histone_ratio_path))
stopifnot(file.exists(Resce_KO_histone_ratio_path))SUZ12 WT proteoforms quantification data.
Code
protein_psi_path <- file.path(tbl_dir_fig4, 'Protein_PSI_3reps_LS.xlsx')3 Main Figure Panels
Import metadata.
Code
mtdt <- read_excel(metadata_path) |>
setNames(c("Sample", "Short_Name", "Clone_Name", "Condition", "Replicate"))
mtdt$Condition <- factor(mtdt$Condition,
levels = c("WT", "CSexon4", "Dexon4", "KO", "KO_L", "KO_S") )3.1 H3 K27-K40 peptide modification ratios
Import DIA data quantification using automatic EpiProfile v2.1 quantifications.
Code
read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) |>
tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set1
read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) |>
# fix an extra _ in the file name
mutate(Sample = gsub(pattern = '2022MH003_NIAR_005_02__30pto_DIA',
replacement = '2022MH003_NIAR_005_02_30pto_DIA',
x = Sample) ) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) |>
tidy_hPTMs(split_double_PTMs = T) -> dia_auto_set2
data_dia_auto <- rbind(dia_auto_set1, dia_auto_set2)Tidy up the information for histone peptide H3 K27- K36. Here I split the H3 K27 and K36 signal coming form the same peptide. In this way I can pile up the all methyl-forms.
Code
data_dia_auto |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
mutate(K27 = case_when( is.na(Sub_PTM1) & grepl(pattern = "^K27", x = Modification) ~ Modification,
grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
mutate(K36 = case_when( is.na(Sub_PTM2) & grepl(pattern = "^K36", x = Modification) ~ Modification,
grepl(pattern = "unmod", x = Modification) ~ "unmod" ) ) |>
relocate(K27, K36, .after = Sub_PTM2) |>
unite(col = "K27", c("Sub_PTM1", "K27"), sep = "", na.rm = T, remove = T ) |>
unite(col = "K36", c("Sub_PTM2", "K36"), sep = "", na.rm = T, remove = T ) |>
mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "unmod",
grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "unmod",
grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
# When one mark is not found I label it unmodified ('unmod') but it could also
# be labelled as "only the other mark".
# mutate(K27 = case_when(grepl(pattern = "^$", x = K27) ~ "Only\nK36me",
# grepl(pattern = "[aA-zZ]", x = K27) ~ K27 )) |>
# mutate(K36 = case_when(grepl(pattern = "^$", x = K36) ~ "Only\nK27ac-me",
# grepl(pattern = "[aA-zZ]", x = K36) ~ K36 )) |>
relocate(K27, K36, .after = PTM) |>
mutate(Histone = factor(Histone, levels = c("H3", "H3.3") ) ) -> data_K27_36Overview of H3 and H3.3 PTMs in lysines 27 and 36
Code
datatable(data_K27_36, rownames = FALSE, filter = 'top',
options = list(pageLength = 7, autoWidth = TRUE) )I can also deconvolve the K27 and K36 signals from the same peptide and pile-up each signal.
Code
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) )|>
select(Histone, Condition, Average_Ratio, K27, Short_Name, Ratio) |>
mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
mutate(Percent = round(Average_Ratio*100, 2) ) |>
select(Condition, Histone, K27, Average_Ratio, Ratio_Sum, Percent) |>
unique() |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
arrange(Condition, Histone, K27) |> ungroup() -> data_K27Same as above but for K36.
Code
data_K27_36 |>
group_by(Condition, K36, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
select(Histone, PTM, Condition, Average_Ratio, K36, Short_Name, Ratio) |>
mutate(Ratio_Sum = sum(Average_Ratio) / 3) |>
mutate(Percent = round(Average_Ratio*100, 2) ) |>
select(Condition, Histone, K36, Average_Ratio, Ratio_Sum, Percent) |>
unique() |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
arrange(Condition, Histone, K36) |> ungroup() -> data_K36For specific histone PTMs one can explore the results in the following tables.
Overview of H3 and H3.3 K27 acetylation and methylation quantifications
Code
datatable(data_K27, rownames = FALSE, filter = 'top',
options = list(pageLength = 5, autoWidth = TRUE) )Overview of H3 and H3.3 K36 methylation quantifications
Code
datatable(data_K36, rownames = FALSE, filter = 'top',
options = list(pageLength = 5, autoWidth = TRUE) )3.1.1 Global H3K27 ratios
Plot figure 4B panel with the pileup of figure H3K27 methylation and acetylation.
Code
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |> # 3 replicates
subset(Histone %in% "H3") |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K27, alluvium = K27,
y = Ratio_Sum, fill = K27) +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
geom_fit_text(aes(label = K27 ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af",
"K27me2" = "#748f46", "K27me3" = "#47632a",
"unmod" = "gray" ),
name = 'H3') +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm -> p_Alluvial_K27
p_Alluvial_K27 
Save to pdf.
Code
ggsave(filename = "Fig4B_H3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4,
plot = p_Alluvial_K27, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)4 Differential abundance analysis of histone PTMs with DEP
Use the area under the peak as counts to be analysed using DEP framework.
Code
rbind( read_EpiProfile_histone_ratio(WT_Dex4_histone_ratio_path),
read_EpiProfile_histone_ratio(Resce_KO_histone_ratio_path) ) |>
# fix an extra _ in the file name like before
mutate(Sample = gsub(pattern = '2022MH003_NIAR_005_02__30pto_DIA',
replacement = '2022MH003_NIAR_005_02_30pto_DIA',
x = Sample) ) |>
left_join(y = mtdt, by = "Sample") |>
mutate(Facility_ID = Sample) |>
mutate(Sample = Short_Name) -> data_dia_auto_not_tidy
Note on the MS peaks areas
One small caveat in this analysis is that the areas under the peaks are very large numbers or at least larger than 2147483647 which is the highest number R can store as integer. To convert the areas to integers from numeric type one could use the package bit64 that stores integers in 64 bits instead of the default 32 bit. My solution was however simpler, I divided the all values by 1000 to scale down by 3 order of magnitude. I then convert the dataframe to matrix.
I tried to add a pseudocounts to the data but it actually doesn’t help with the analysis and is better to input the areas as they are. The function make_unique is the first step of the analysis of DEP and creates a unique data frame.
Code
data_dia_auto_not_tidy |>
tidy_hPTMs() |>
select(First_Col, PTM, Sample, Area) |>
# reduce all areas by a 1000 to avoid having to coerc to integer big numbers
mutate(Area = Area / 1e4) |>
pivot_wider(id_cols = c(First_Col, PTM), names_from = Sample, values_from = Area) |>
# Coerce the areas to integers
mutate(across(.cols = ends_with(c("_1", "_2", "_3")), .fns = as.integer )) |>
make_unique(names = "First_Col", ids = "PTM") -> data_uniqueThe data frame data_unique contains 22 columns (of which 18 are samples) and 185 rows corresponding to the identified histone PTMs.
Create a metadata compatible for DEP set the KO as the reference sample as I’m interesting in the rescue efficiency relative to the KO.
Code
# Set metadata colnames to be as `make_se` wants them to be
colnames(mtdt) <- c("Sample", "label", "Clone_Name", "condition", "replicate")
mutate(mtdt, condition = factor(condition, levels = c("WT", "CSexon4", "Dexon4",
"KO", "KO_L", "KO_S"))) |>
print(n = 3) -> mtdt_dep# A tibble: 18 × 5
Sample label Clone_Name condition replicate
<chr> <chr> <chr> <fct> <dbl>
1 2022MZ006_IVMO_001_02_25pto_DIA WT_1 ZWT WT 1
2 2022MZ006_IVMO_002_02_25pto_DIA WT_2 ZA1 WT 2
3 2022MZ006_IVMO_003_02_25pto_DIA WT_3 ZA2 WT 3
# ℹ 15 more rowsFind the columns that contain the samples areas in the data_unique dataframe and create a SummarizedExperiment object to be used in for the DEP analysis. The first step is to do a log2 transformation of the area counts.
Code
Samples_area_columns <- grep(pattern = paste(mtdt_dep$label, collapse = "|"), colnames(data_unique))
data_se <- make_se(data_unique, columns = Samples_area_columns, expdesign = mtdt_dep)Now the data_se contains the modifications in the rows (assuming they are proteins) and in the column the area.
Code
data_filt <- filter_missval(data_se, thr = 0)Check PTMs frequency across samples.
Code
plot_frequency(data_filt)
Overall most modifications are found in all samples.
Code
data_norm <- normalize_vsn(data_filt)Verify the fit.
Code
meanSdPlot(x = data_norm)
Seems okay.
Check normalisation effect on data.
Code
plot_normalization(data_se, data_filt, data_norm) +
scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
"CSexon4" = "#B0461C", "KO" = "gray45",
"KO_L" = "mediumpurple3", "KO_S" = "darkorange1")) +
theme(axis.text = element_text(color = 'black'),
plot.background = element_blank(),
panel.background = element_blank())
Check missing values.
Code
plot_missval(data_norm)
I don’t see obvious missing values trends, besides only that replicate #3 of CSexon has less PTMs detected. I don’t think it’s necessary to impute missing values.
4.1 Test for significance against WT control
Test every sample versus WT control
Code
data_diff <- test_diff(data_norm, type = "control", control = "WT")Significance thresholds:
P-value <= 0.05
|Log2 Fold Change| >= 1.5
Code
dep <- add_rejections(data_diff, alpha = 0.05, lfc = 1.5 )Check sample clustering a Gower’s distance matrix.s
Code
plot_dist(dep, significant = TRUE, pal = "PuOr") 
This is basically showing the 2 batches of the 2 MS acquisition experiments.
Check sample clustering using the same distance.
Code
plot_heatmap(dep, type = "centered", kmeans = F, clustering_distance = "gower",
col_limit = 4, show_row_names = TRUE,
indicate = c("condition", "replicate"))
Get results to make a volcano plot.
Code
data_results <- get_results(dep)
data_results |>
as_tibble() |>
select(-significant) |>
pivot_longer(cols = !c("name", "ID"),
names_to = c("contrast", "Feature"),
names_sep = c("_(vs_WT_|centered)"),
values_to = "Value"
) |>
mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
pivot_wider(names_from = "Feature", values_from = "Value") |>
dplyr::select(-ID) -> res_all_tidy
colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'Extract WT sample area of each PTM, to be used for the volcano plot.
Code
data_dia_auto |>
subset(Condition == "WT") |>
group_by(First_Col) |>
mutate(WT_Mean_Area = mean(Area) + 1 ) |>
relocate(WT_Mean_Area, .after = PTM) |>
select(c(First_Col, WT_Mean_Area)) |>
unique() -> WT_areaGet the histone PTMs info columns.
Code
histone_ptms_anno <- unique(data_dia_auto[, 1:10])Add all the info the the tidy results dataframe.
Code
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )
res_all |>
mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
arrange(desc(Histone), Peptide_Start, Peptide_End) |>
ungroup() -> res_allExport all tidy results for supplementary table 3 as tab delimited, or as csv with UTF-8 Byte order mark which indicates to Excel the csv is UTF-8 encoded.
Code
write_delim(x = res_all, file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.tab'),
delim = '\t', col_names = T, append = F, quote = 'none', na = "NA",
progress = F, escape = 'none')
write_excel_csv(x = res_all,
file = file.path(tbl_dir_fig4, 'DEP_analysis_histone_MS.csv'),
delim = ',', col_names = T, append = F, quote = 'none',
na = "NA", progress = F, escape = 'none')Interactively explore the table
Code
res_all |>
mutate(across(.cols = where(is.numeric), \(x) round(x, 4) ) ) |>
datatable(rownames = FALSE, filter = 'top',
options = list(pageLength = 3, autoWidth = TRUE) )4.2 PTMs centred abundances
Plot centred abundance
Code
res_all |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
subset(!is.na(centered)) |>
ggplot(aes(x = centered, y = PTM, fill = contrast)) +
facet_grid(rows = vars(Modification), scales = 'free' ) +
geom_vline(xintercept = 0, linewidth = 0.4) +
geom_col(width = 0.75,
position = position_dodge(width = 0.75, preserve = 'single')) +
geom_point(aes(size = WT_Mean_Area),
shape = 21, position = position_dodge(width = 0.75),
stroke = 0.2) +
labs(x = "Centered PTM abundance") +
scale_x_continuous(limits = c(-5, 5), oob = scales::oob_squish ) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
scale_fill_manual(values = c("WT" = "royalblue3", "Dexon4" = "goldenrod",
"CSexon4" = "#B0461C", "KO" = "gray45",
"KO_L" = "mediumpurple3", "KO_S" = "darkorange1"),
name = "Sample") +
theme_classic() +
theme(strip.text = element_blank(),
strip.background = element_blank(),
axis.text = element_text(colour = "black"),
axis.title.y = element_blank(),
legend.key.size = unit(1, 'mm'),
panel.grid.major.y = element_line()) -> p_centered_K27_K36
p_centered_K27_K36
4.3 Volcano plot
Show differential histone PTMs as volcano plots.
4.3.1 ∆ex4 vs WT
Prepare data for plot.
Code
res_Dex4 <- subset(res_all, contrast == "Dexon4" & !is.na(p.val))
res_Dex4 <- res_Dex4 |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot figure 4C.
Code
ggplot(res_Dex4 ) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 5.5, y = 0.10, label = "∆ex4",
colour = "black", fill = "goldenrod", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -5.5, y = 0.10, label = "WT",
colour = "white", fill = "royalblue3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_Dex4, p.adj < 0.09 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_Dex4, p.adj < 0.1 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-6.4, 6.4), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.91, 0.87),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_Dex4
p_Dex4
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "Fig4C_Dex4_vs_WT_Volcano_w_labels.pdf",
plot = p_Dex4, device = cairo_pdf, units = "cm",
height = 5.4, width = 7.0)4.3.2 CSex4 vs WT
Prepare data of CSex4 vs WT
Code
res_CSex4 <- subset(res_all, contrast == "CSexon4" & !is.na(p.val))
res_CSex4 <- res_CSex4 |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Note: nothing is differentially expressed in this contrast
Code
table(res_CSex4$Direction)
None
159 4.3.3 Volcano plot KO vs WT
Code
res_KO <- subset(res_all, contrast == "KO" & !is.na(p.val))
res_KO <- res_KO |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) Plot KO vs WT
Code
ggplot(res_KO ) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 8.5, y = 0.10, label = "KO",
colour = "white", fill = "gray16", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -12.5, y = 0.10, label = "WT",
colour = "white", fill = "royalblue3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_KO, p.adj <= 0.15 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-15, 10), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.25), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.91, 0.67),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_KO
p_KO
Note: these WT and KO samples were processed at different moment, so there could be some batch effect difference that limits the statistical power.
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_KO_vs_WT_Volcano_w_labels.pdf",
plot = p_KO, device = cairo_pdf, units = "cm",
height = 12, width = 10)4.4 Test for significance vs the KO+L
Test against Long Suz12 isoform KO rescue.
Code
data_diff_long <- test_diff(data_norm, type = "control", control = "KO_L")Use same significance thresholds as before.
Code
dep_rescues <- add_rejections(data_diff_long, alpha = 0.05, lfc = 1.5 )Get results to make a volcano plot of the long vs the short
Code
data_results <- get_results(dep_rescues)
data_results |>
as_tibble() |>
select(-significant) |>
pivot_longer(cols = !c("name", "ID"),
names_to = c("contrast", "Feature"),
names_sep = c("_(vs_KO_L_|centered)"),
values_to = "Value"
) |>
mutate(Feature = ifelse(test = Feature == "", yes = "centered", no = Feature)) |>
pivot_wider(names_from = "Feature", values_from = "Value") |>
dplyr::select(-ID) -> res_all_tidy
colnames(res_all_tidy)[colnames(res_all_tidy) == "name"] <- 'First_Col'Add all the info the the tidy results dataframe.
Code
res_all <- left_join(x = res_all_tidy, y = histone_ptms_anno, by = join_by(First_Col))
res_all <- left_join(res_all, WT_area, by = join_by(First_Col) )
res_all |>
mutate(Histone = factor(Histone, levels = sort(unique(res_all$Histone)) ) ) |>
arrange(desc(Histone), Peptide_Start, Peptide_End) |>
ungroup() -> res_allPrepare data for plot.
Code
res_KOrS <- subset(res_all, contrast == "KO_S" & !is.na(p.val))
res_KOrS <- res_KOrS |>
mutate(Direction = case_when(p.adj <= 0.05 & ratio > 0 ~ 'UP',
p.adj <= 0.05 & ratio < 0 ~ 'DOWN',
p.adj > 0.05 ~ 'None') ) 4.4.1 Volcano plot rescue SUZ12-L vs SUZ12-S
Plot rescues.
Code
ggplot(res_KOrS) +
aes(x = ratio, y = -log10(p.val), fill = Direction, size = WT_Mean_Area ) +
geom_point(shape = 21, stroke = 0.2) +
annotate(geom = "label", x = 7.5, y = 0.10, label = "KO+S",
colour = "black", fill = "darkorange1", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
annotate(geom = "label", x = -4.85, y = 0.10, label = "KO+L",
colour = "white", fill = "mediumpurple3", size = 2,
label.padding = grid::unit(0.5, "mm"),
label.r = unit(0.25, "mm"), family = "Arial",
label.size = grid::unit(0.125, "mm") ) +
# points on the right side
geom_label_repel(data = subset(res_KOrS, p.adj < 0.1 & ratio > 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = -0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
# points on the left side
geom_label_repel(data = subset(res_KOrS, p.adj < 0.1 & ratio < 0 ),
aes(label = PTM), fill = 'white',
seed = 16, show.legend = F, segment.curvature = -1e-20,
family = "Arial", size = 1.85, nudge_x = 0.25,
segment.color = 'black',verbose = F,
box.padding = grid::unit(1, "mm"),
point.padding = grid::unit(0.55, "mm"),
label.padding = grid::unit(0.5, "mm"),
label.size = 0.120,
max.overlaps = 20) +
labs(x = expression(log[2] ~ "Fold Change"), y = expression(-log[10] ~ "P-Value")) +
scale_fill_manual(values = c('dodgerblue', "gray84", "firebrick3"), guide = 'none') +
guides(size = guide_legend(override.aes = list(fill = "white"))) +
scale_x_continuous(expand = expansion(add = 0.02, mult = 0), limits = c(-5.5, 8.25), n.breaks = 7) +
scale_y_continuous(expand = expansion(add = c(0, 0.05), mult = 0), limits = c(0, NA)) +
scale_size_area(breaks = c(1e8, 1e9, 1e10, 1e11, 1e12),
labels = scientific_10, name = "WT PTM\nabundance") +
theme_classic(base_size = 6, base_family = "Arial") +
theme(panel.background = element_blank(),
panel.border = element_blank(),
panel.grid.major = element_line(linewidth = 0.2),
legend.position = c(0.10, 0.83),
legend.key = element_blank(),
legend.key.size = unit(1.0, units = 'mm'),
legend.text = element_text(size = 5, margin = margin(l = -1, r = -1, unit = 'mm')),
legend.direction = "vertical",
legend.background = element_blank(),
legend.box.background = element_blank(),
legend.box.spacing = unit(-3, 'mm'),
axis.text = element_text(colour = 'black'),
axis.ticks.length = unit(1, units = 'mm'),
axis.ticks = element_line(colour = 'black'),
axis.title = element_text(size = 5),
axis.line = element_line(linewidth = 0.2),
plot.background = element_blank() ) -> p_rescues
p_rescues
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_KOrS_vs_KOrL_Volcano_w_labels.pdf",
plot = p_rescues, device = cairo_pdf, units = "cm",
height = 6, width = 7.9)
Things I could still add to this analysis
I could still add:
rescue efficiency
heatmap with all PTMs.
Now plot to the supplementary figures
5 Supplementary Figure Panels
5.1 K27-K40 peptide modification ratios
5.1.1 H3.3K27
There’s no H3.3 signal in the main figure, here I plot only H3.3 signal.
Code
data_K27_36 |>
group_by(Condition, K27, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
subset(Histone %in% "H3.3") |>
mutate(K27 = factor(K27, levels = c("K27ac", "K27me3", "K27me2", "K27me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K27, alluvium = K27,
y = Ratio_Sum, fill = K27) +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, colour = "black") +
geom_fit_text(aes(label = K27 ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = c("K27ac" = "#feac81", "K27me1" ="#ced1af",
"K27me2" = "#748f46", "K27me3" = "#47632a",
"unmod" = "gray" ),
name = 'H3.3') +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm -> p_Alluvial_H33K27
p_Alluvial_H33K27 
Save to pdf figure S4D.
Code
ggsave(filename = "FigS4D_H3.3_only_K27_Mrgd_Alluvial.pdf", path = pdf_dir_fig4,
plot = p_Alluvial_H33K27, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)5.1.2 H3 and H3.3 K36
Plot the pileup of H3 & H3.3 K36 methylation.
Code
data_K27_36 |>
group_by(Condition, K36, Histone) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3) |>
subset(Histone == c("H3", "H3.3") ) |>
mutate(K36 = factor(K36, levels = c("K36me3", "K36me2", "K36me1", "unmod") ) ) |>
ggplot() +
aes(x = Condition, stratum = K36, alluvium = K36,
y = Ratio_Sum, fill = K36) +
facet_wrap(~ Histone, scales = 'free_y') +
geom_alluvium(alpha = 0.5, width = 0.4) +
geom_stratum(alpha = 1, width = 0.6, linewidth = 0.1, colour = "black") +
geom_fit_text(aes(label = K36 ), stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"), size = 7, family = "Arial",
padding.y = grid::unit(0.2, "mm"),
show.legend = F) +
scale_fill_manual(values = c("K36me1" = "#87B7AA", "K36me2" = "#8FA9C2",
"K36me3" = "#8386B7", "unmod" = "gray" ) ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO', 'KO+L', 'KO+S')
) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.title = element_blank()) -> p_Alluvial_K36
p_Alluvial_K36
Save to pdf figure S4E
Code
ggsave(path = pdf_dir_fig4, filename = "FigS4E_H3_H3.3_K36_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K36, device = cairo_pdf, units = "cm",
width = 12.0, height = 5)5.2 Extra plots
Here I plot extra histone PTMs that were not included in the publication in case someone is interested.
5.2.1 H3K4
Plot H3K4 modifications.
Code
data_dia_auto |>
group_by(Condition, Modification) |>
subset(Peptide_Start >= 3 & Peptide_End <= 8) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Demuth", direction = 1, n = 5), name = 'H3' ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
guides(fill = guide_legend(byrow = TRUE) ) + # required for the legend item spacing
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -9, unit = "mm")),
legend.key.width = unit(6.5, "mm") ) -> p_Alluvial_K4
p_Alluvial_K4
Save to pdf.
Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3K4_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K4, device = cairo_pdf, units = "cm",
width = 6.5, height = 4.85)5.2.2 H3K9
Plot all combinatorial K9 and K14 modifications on histone H3.
Code
data_dia_auto |>
group_by(Condition, Modification) |>
subset(Histone == "H3") |>
subset(Peptide_Start >= 9 & Peptide_End <= 17) |>
mutate(Average_Ratio = mean(Ratio) ) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum2(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Hokusai1", direction = -1, n = 15), name = 'H3' ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S')) +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -12.0, unit = "mm")),
legend.key.width = unit(11.2, "mm"),
legend.key.height = unit(1, 'mm'),
legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K9
p_Alluvial_K9
Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3K9_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K9, device = cairo_pdf, units = "cm",
width = 6.4, height = 4.85)5.2.3 All combinations of H3/H3.3 K27 & K36
Plot all combinatorial K27 and K36 modifications on histone H3 and H3.3 variant.
Code
data_dia_auto |>
group_by(Condition, Modification, Histone) |>
subset(Peptide_Start >= 27 & Peptide_End <= 40) |>
mutate(Average_Ratio = mean(Ratio)) |>
summarise(Ratio_Sum = sum(Average_Ratio) / 3, .groups = 'keep') |>
ggplot() +
aes(x = Condition, stratum = Modification, alluvium = Modification,
y = Ratio_Sum, fill = Modification) +
facet_wrap(~ Histone) +
geom_alluvium(alpha = 0.4, width = 0.4) +
geom_stratum(alpha = 1, width = 0.6, size = 0.1) +
geom_fit_text(aes(label = Modification ),
stat = "stratum", width = 0.6, min.size = 2,
padding.x = grid::unit(0.2, "mm"),
padding.y = grid::unit(0.2, "mm"),
size = 7, family = "Arial", show.legend = F) +
scale_fill_manual(values = met.brewer(name = "Hiroshige", direction = 1, n = 15) ) +
scale_x_discrete(expand = expansion(mult = 0.05, add = c(0, 0.1)),
labels = c('WT', 'CSex4', '∆ex4', 'KO',
'KO+L', 'KO+S'), name = 'H3 & H3.3') +
scale_y_continuous(expand = expansion(mult = 0, add = 0) ) +
labs(y = "Average Ratio") +
alluvial_thm +
theme(legend.text = element_text(margin = margin(l = -14.5, unit = "mm")),
legend.key.width = unit(13.75, "mm"),
legend.key.height = unit(1, 'mm'),
legend.spacing.x = unit(1, 'mm') ) -> p_Alluvial_K27_K36
p_Alluvial_K27_K36
This is the best representative picture of K27-K40 modifications relative abundance. This also confirms how the H3.3K36me3 is really the most dramatic change between KO+L and KO+S.
Code
ggsave(path = pdf_dir_fig4, filename = "FigExtra_H3_H3.3_K27_K36_Mrgd_Alluvial.pdf",
plot = p_Alluvial_K27_K36, device = cairo_pdf, units = "cm",
width = 14, height = 5)End analysis.
5.3 Protein PSI quantification
Code
Protein_PSI_3reps_LS <- read_excel(path = protein_psi_path, sheet = "Sheet1") |>
pivot_longer(cols = starts_with('SUZ12'), names_to = 'Proteoform', values_to = 'PSI') |>
mutate(PSI = PSI * 100)
Protein_PSI_3reps_LS |> group_by(Proteoform) |> summarise(median(PSI))# A tibble: 2 × 2
Proteoform `median(PSI)`
<chr> <dbl>
1 SUZ12-L 90.6
2 SUZ12-S 9.43Plot
Code
ggplot(Protein_PSI_3reps_LS) +
aes(x = Proteoform, y = PSI, fill = Proteoform) +
geom_boxplot(show.legend = F, outlier.shape = NA, width = 0.5, lwd = 0.1) +
geom_point(shape = 21, size = 0.85, stroke = 0.1, show.legend = F) +
coord_cartesian(ylim = c(0, 100), clip = 'off') +
scale_fill_manual(values = c("SUZ12-L" = "mediumpurple3", "SUZ12-S" = "darkorange1")) +
scale_y_continuous(expand = expansion(mult = 0, add = 0),
n.breaks = 10) +
labs(x = 'Proteoform',
y = 'Rel. abundance in WT ESCs') +
alluvial_thm +
theme(panel.grid.major.y = element_line(colour = 'gray84', linewidth = 0.2),
axis.title.x = element_text(size = 5) ) -> p_prot_PSI
p_prot_PSI
Save to pdf figure S4C
Code
ggsave(filename = "FigS4C_SUZ12_WT_proteoforms_quantification_barplot.pdf",
plot = p_prot_PSI, device = cairo_pdf, path = pdf_dir_fig4, units = 'cm',
width = 2.65, height = 2.75)6 Session Info
Code
sessioninfo::session_info()─ Session info ───────────────────────────────────────────────────────────────
setting value
version R version 4.1.2 (2021-11-01)
os macOS Catalina 10.15.7
system x86_64, darwin17.0
ui X11
language (EN)
collate C
ctype UTF-8
tz Europe/Madrid
date 2024-02-07
pandoc 2.10.1 @ /usr/local/bin/ (via rmarkdown)
─ Packages ───────────────────────────────────────────────────────────────────
! package * version date (UTC) lib source
ade4 1.7-22 2023-02-06 [1] CRAN (R 4.1.2)
affy 1.72.0 2021-10-26 [1] Bioconductor
affyio 1.64.0 2021-10-26 [1] Bioconductor
annotate 1.72.0 2021-10-26 [1] Bioconductor
AnnotationDbi 1.56.2 2021-11-09 [1] Bioconductor
assertthat 0.2.1 2019-03-21 [1] CRAN (R 4.1.0)
Biobase 2.54.0 2021-10-26 [1] Bioconductor
BiocFileCache 2.2.1 2022-01-23 [1] Bioconductor
BiocGenerics 0.40.0 2021-10-26 [1] Bioconductor
BiocManager 1.30.22 2023-08-08 [1] CRAN (R 4.1.2)
BiocParallel 1.28.3 2021-12-09 [1] Bioconductor
biomaRt 2.50.3 2022-02-06 [1] Bioconductor
Biostrings 2.62.0 2021-10-26 [1] Bioconductor
bit 4.0.5 2022-11-15 [1] CRAN (R 4.1.2)
bit64 4.0.5 2020-08-30 [1] CRAN (R 4.1.0)
bitops 1.0-7 2021-04-24 [1] CRAN (R 4.1.0)
blob 1.2.4 2023-03-17 [1] CRAN (R 4.1.2)
bslib 0.5.1 2023-08-11 [1] CRAN (R 4.1.2)
cachem 1.0.8 2023-05-01 [1] CRAN (R 4.1.2)
Cairo * 1.6-1 2023-08-18 [1] CRAN (R 4.1.2)
callr 3.7.3 2022-11-02 [1] CRAN (R 4.1.2)
cellranger 1.1.0 2016-07-27 [1] CRAN (R 4.1.0)
circlize 0.4.15 2022-05-10 [1] CRAN (R 4.1.2)
cli 3.6.1 2023-03-23 [1] CRAN (R 4.1.2)
clue 0.3-64 2023-01-31 [1] CRAN (R 4.1.2)
cluster 2.1.4 2022-08-22 [1] CRAN (R 4.1.2)
codetools 0.2-19 2023-02-01 [1] CRAN (R 4.1.2)
colorspace 2.1-0 2023-01-23 [1] CRAN (R 4.1.2)
ComplexHeatmap 2.10.0 2021-10-26 [1] Bioconductor
crayon 1.5.2 2022-09-29 [1] CRAN (R 4.1.2)
crosstalk 1.2.0 2021-11-04 [1] CRAN (R 4.1.0)
csaw 1.28.0 2021-10-26 [1] Bioconductor
curl 5.2.0 2023-12-08 [1] CRAN (R 4.1.2)
DBI 1.1.3 2022-06-18 [1] CRAN (R 4.1.2)
dbplyr 2.3.3 2023-07-07 [1] CRAN (R 4.1.2)
DelayedArray 0.20.0 2021-10-26 [1] Bioconductor
DEP * 1.16.0 2021-10-26 [1] Bioconductor
desc 1.4.2 2022-09-08 [1] CRAN (R 4.1.2)
DESeq2 1.34.0 2021-10-26 [1] Bioconductor
devtools 2.4.5 2022-10-11 [1] CRAN (R 4.1.2)
digest 0.6.33 2023-07-07 [1] CRAN (R 4.1.2)
doParallel 1.0.17 2022-02-07 [1] CRAN (R 4.1.2)
dplyr * 1.1.2 2023-04-20 [1] CRAN (R 4.1.2)
DT * 0.29 2023-08-29 [1] CRAN (R 4.1.2)
edgeR 3.36.0 2021-10-26 [1] Bioconductor
ellipsis 0.3.2 2021-04-29 [1] CRAN (R 4.1.0)
evaluate 0.21 2023-05-05 [1] CRAN (R 4.1.2)
fansi 1.0.4 2023-01-22 [1] CRAN (R 4.1.2)
farver 2.1.1 2022-07-06 [1] CRAN (R 4.1.2)
fastmap 1.1.1 2023-02-24 [1] CRAN (R 4.1.2)
fdrtool 1.2.17 2021-11-13 [1] CRAN (R 4.1.0)
filelock 1.0.2 2018-10-05 [1] CRAN (R 4.1.0)
forcats 1.0.0 2023-01-29 [1] CRAN (R 4.1.2)
foreach 1.5.2 2022-02-02 [1] CRAN (R 4.1.2)
foreign 0.8-84 2022-12-06 [1] CRAN (R 4.1.2)
fs 1.6.3 2023-07-20 [1] CRAN (R 4.1.2)
genefilter 1.76.0 2021-10-26 [1] Bioconductor
geneplotter 1.72.0 2021-10-26 [1] Bioconductor
generics 0.1.3 2022-07-05 [1] CRAN (R 4.1.2)
GenomeInfoDb 1.30.1 2022-01-30 [1] Bioconductor
GenomeInfoDbData 1.2.7 2023-08-29 [1] Bioconductor
GenomicRanges 1.46.1 2021-11-18 [1] Bioconductor
GetoptLong 1.0.5 2020-12-15 [1] CRAN (R 4.1.0)
ggalluvial 0.12.5 2023-02-22 [1] CRAN (R 4.1.2)
ggfittext 0.10.0 2023-04-04 [1] CRAN (R 4.1.2)
ggplot2 3.4.3 2023-08-14 [1] CRAN (R 4.1.2)
ggrepel 0.9.3 2023-02-03 [1] CRAN (R 4.1.2)
ggseqlogo 0.1 2017-07-25 [1] CRAN (R 4.1.0)
GlobalOptions 0.1.2 2020-06-10 [1] CRAN (R 4.1.0)
glue 1.6.2 2022-02-24 [1] CRAN (R 4.1.2)
gmm 1.8 2023-06-06 [1] CRAN (R 4.1.2)
gtable 0.3.4 2023-08-21 [1] CRAN (R 4.1.2)
hexbin 1.28.3 2023-03-21 [1] CRAN (R 4.1.2)
hms 1.1.3 2023-03-21 [1] CRAN (R 4.1.2)
htmltools 0.5.6 2023-08-10 [1] CRAN (R 4.1.2)
htmlwidgets 1.6.2 2023-03-17 [1] CRAN (R 4.1.2)
httpuv 1.6.11 2023-05-11 [1] CRAN (R 4.1.2)
httr 1.4.7 2023-08-15 [1] CRAN (R 4.1.2)
impute 1.68.0 2021-10-26 [1] Bioconductor
imputeLCMD 2.1 2022-06-10 [1] CRAN (R 4.1.2)
IRanges 2.28.0 2021-10-26 [1] Bioconductor
iterators 1.0.14 2022-02-05 [1] CRAN (R 4.1.2)
jquerylib 0.1.4 2021-04-26 [1] CRAN (R 4.1.0)
jsonlite 1.8.7 2023-06-29 [1] CRAN (R 4.1.2)
KEGGREST 1.34.0 2021-10-26 [1] Bioconductor
knitr 1.43 2023-05-25 [1] CRAN (R 4.1.2)
labeling 0.4.2 2020-10-20 [1] CRAN (R 4.1.0)
later 1.3.1 2023-05-02 [1] CRAN (R 4.1.2)
lattice 0.21-8 2023-04-05 [1] CRAN (R 4.1.2)
lifecycle 1.0.3 2022-10-07 [1] CRAN (R 4.1.2)
limma 3.50.3 2022-04-07 [1] Bioconductor
locfit 1.5-9.8 2023-06-11 [1] CRAN (R 4.1.2)
magick 2.7.5 2023-08-07 [1] CRAN (R 4.1.2)
magrittr 2.0.3 2022-03-30 [1] CRAN (R 4.1.2)
MALDIquant 1.22.1 2023-03-20 [1] CRAN (R 4.1.2)
MASS 7.3-60 2023-05-04 [1] CRAN (R 4.1.2)
Matrix 1.6-1 2023-08-14 [1] CRAN (R 4.1.2)
MatrixGenerics 1.6.0 2021-10-26 [1] Bioconductor
matrixStats 1.0.0 2023-06-02 [1] CRAN (R 4.1.2)
memoise 2.0.1 2021-11-26 [1] CRAN (R 4.1.0)
metapod 1.2.0 2021-10-26 [1] Bioconductor
MetBrewer 0.2.0 2022-03-21 [1] CRAN (R 4.1.2)
mime 0.12 2021-09-28 [1] CRAN (R 4.1.0)
miniUI 0.1.1.1 2018-05-18 [1] CRAN (R 4.1.0)
msa 1.26.0 2021-10-26 [1] Bioconductor
MsCoreUtils 1.6.2 2022-02-24 [1] Bioconductor
MSnbase 2.20.4 2022-01-16 [1] Bioconductor
munsell 0.5.0 2018-06-12 [1] CRAN (R 4.1.0)
mvtnorm 1.2-3 2023-08-25 [1] CRAN (R 4.1.2)
mzID 1.32.0 2021-10-26 [1] Bioconductor
mzR 2.28.0 2021-10-27 [1] Bioconductor
ncdf4 1.21 2023-01-07 [1] CRAN (R 4.1.2)
R niar * 0.3.0 <NA> [?] <NA>
norm 1.0-11.1 2023-06-18 [1] CRAN (R 4.1.2)
patchwork 1.1.3 2023-08-14 [1] CRAN (R 4.1.2)
pcaMethods 1.86.0 2021-10-26 [1] Bioconductor
pheatmap * 1.0.12 2019-01-04 [1] CRAN (R 4.1.0)
pillar 1.9.0 2023-03-22 [1] CRAN (R 4.1.2)
pkgbuild 1.4.2 2023-06-26 [1] CRAN (R 4.1.2)
pkgconfig 2.0.3 2019-09-22 [1] CRAN (R 4.1.0)
pkgload 1.3.2.1 2023-07-08 [1] CRAN (R 4.1.2)
plyr 1.8.8 2022-11-11 [1] CRAN (R 4.1.2)
png 0.1-8 2022-11-29 [1] CRAN (R 4.1.2)
preprocessCore 1.56.0 2021-10-26 [1] Bioconductor
prettyunits 1.1.1 2020-01-24 [1] CRAN (R 4.1.0)
processx 3.8.2 2023-06-30 [1] CRAN (R 4.1.2)
profvis 0.3.8 2023-05-02 [1] CRAN (R 4.1.2)
progress 1.2.2 2019-05-16 [1] CRAN (R 4.1.0)
promises 1.2.1 2023-08-10 [1] CRAN (R 4.1.2)
ProtGenerics 1.26.0 2021-10-26 [1] Bioconductor
ps 1.7.5 2023-04-18 [1] CRAN (R 4.1.2)
psychTools 2.3.6 2023-06-18 [1] CRAN (R 4.1.2)
purrr 1.0.2 2023-08-10 [1] CRAN (R 4.1.2)
R6 2.5.1 2021-08-19 [1] CRAN (R 4.1.0)
rappdirs 0.3.3 2021-01-31 [1] CRAN (R 4.1.0)
RColorBrewer 1.1-3 2022-04-03 [1] CRAN (R 4.1.2)
Rcpp 1.0.11 2023-07-06 [1] CRAN (R 4.1.2)
RCurl 1.98-1.12 2023-03-27 [1] CRAN (R 4.1.2)
readr * 2.1.4 2023-02-10 [1] CRAN (R 4.1.2)
readxl * 1.4.3 2023-07-06 [1] CRAN (R 4.1.2)
remotes 2.4.2.1 2023-07-18 [1] CRAN (R 4.1.2)
rjson 0.2.21 2022-01-09 [1] CRAN (R 4.1.2)
rlang 1.1.1 2023-04-28 [1] CRAN (R 4.1.2)
rmarkdown 2.24 2023-08-14 [1] CRAN (R 4.1.2)
rprojroot 2.0.3 2022-04-02 [1] CRAN (R 4.1.2)
Rsamtools 2.10.0 2021-10-26 [1] Bioconductor
RSQLite 2.3.1 2023-04-03 [1] CRAN (R 4.1.2)
rstudioapi 0.15.0 2023-07-07 [1] CRAN (R 4.1.2)
S4Vectors 0.32.4 2022-03-29 [1] Bioconductor
sandwich 3.0-2 2022-06-15 [1] CRAN (R 4.1.2)
sass 0.4.7 2023-07-15 [1] CRAN (R 4.1.2)
scales 1.2.1 2022-08-20 [1] CRAN (R 4.1.2)
seqinr 4.2-30 2023-04-05 [1] CRAN (R 4.1.2)
sessioninfo 1.2.2 2021-12-06 [1] CRAN (R 4.1.0)
shape 1.4.6 2021-05-19 [1] CRAN (R 4.1.0)
shiny 1.7.5 2023-08-12 [1] CRAN (R 4.1.2)
shinydashboard 0.7.2 2021-09-30 [1] CRAN (R 4.1.0)
stringi 1.7.12 2023-01-11 [1] CRAN (R 4.1.2)
stringr 1.5.0 2022-12-02 [1] CRAN (R 4.1.2)
SummarizedExperiment 1.24.0 2021-10-26 [1] Bioconductor
survival 3.5-7 2023-08-14 [1] CRAN (R 4.1.2)
tibble 3.2.1 2023-03-20 [1] CRAN (R 4.1.2)
tidyr 1.3.0 2023-01-24 [1] CRAN (R 4.1.2)
tidyselect 1.2.0 2022-10-10 [1] CRAN (R 4.1.2)
tmvtnorm 1.5 2022-03-22 [1] CRAN (R 4.1.2)
tzdb 0.4.0 2023-05-12 [1] CRAN (R 4.1.2)
urlchecker 1.0.1 2021-11-30 [1] CRAN (R 4.1.0)
usethis 2.2.2 2023-07-06 [1] CRAN (R 4.1.2)
utf8 1.2.3 2023-01-31 [1] CRAN (R 4.1.2)
vctrs 0.6.3 2023-06-14 [1] CRAN (R 4.1.2)
vroom 1.6.3 2023-04-28 [1] CRAN (R 4.1.2)
vsn 3.62.0 2021-10-26 [1] Bioconductor
withr 2.5.0 2022-03-03 [1] CRAN (R 4.1.2)
xfun 0.40 2023-08-09 [1] CRAN (R 4.1.2)
XICOR 0.4.1 2023-04-21 [1] CRAN (R 4.1.2)
XML 3.99-0.14 2023-03-19 [1] CRAN (R 4.1.2)
xml2 1.3.5 2023-07-06 [1] CRAN (R 4.1.2)
xtable 1.8-4 2019-04-21 [1] CRAN (R 4.1.0)
XVector 0.34.0 2021-10-26 [1] Bioconductor
yaml 2.3.7 2023-01-23 [1] CRAN (R 4.1.2)
zlibbioc 1.40.0 2021-10-26 [1] Bioconductor
zoo 1.8-12 2023-04-13 [1] CRAN (R 4.1.2)
[1] /Library/Frameworks/R.framework/Versions/4.1/Resources/library
R ── Package was removed from disk.
──────────────────────────────────────────────────────────────────────────────References
Yuan, Zuo-Fei, Simone Sidoli, Dylan M Marchione, Johayra Simithy, Kevin A Janssen, Mary R Szurgot, and Benjamin A Garcia. 2018. “EpiProfile 2.0: A Computational Platform for Processing Epi-Proteomics Mass Spectrometry Data.” Journal of Proteome Research 17 (7): 2533–41. https://doi.org/10.1021/acs.jproteome.8b00133.
Zhang, Xiaofei, Arne H. Smits, Gabrielle BA van Tilburg, Huib Ovaa, Wolfgang Huber, and Michiel Vermeulen. 2018. “Proteome-wide identification of ubiquitin interactions using UbIA-MS.” Nature Protocols 13 (3): 530–50. https://doi.org/10.1038/nprot.2017.147.